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1 Introduction 

The computation of the Delaunay triangulation of a set of n points in the 
plane is one of the classical problems in computational geometry and plenty of 
algorithms have been proposed to solve it. 

These Delaunay algorithms can have different characteristics: 

• Optimal on worst case data, i.e. O(nlogn) time. 

• Good complexity on random data only 

• Randomized 

• On-line vs off-line 

In the current trade-off between algorithmic simplicity, practical efficiency 
and theoretical optimality, practitioners often choose the simplicity and practical 
efficiency taking the risk of having bad performance on some special kind of data. 

Our aim is to conciliate many of the above aspects, namely to obtain an 
incremental algorithm using simple data structure having good practical per- 
formance on realistic input and still provable O(nlogn) computation time on 
any data set. 

Previous related work 

Our work is strongly related to some previous algorithms for Delaunay trian- 
gulation. All these algorithms are incremental and their complexity is random- 
ized, they use some location structure to find where the new point is inserted, 
and then update the triangulation. 

The first idea of a randomized incremental construction for the Delaunay 



triangulation [BT86| uses a location structure based on the history of the De- 
launay triangulation: the Delaunay tree. Point Pi is inserted at time i, and to 
find where point p n fell, p n is located in all the triangulations at times 1 to n— 1; 
the location at time i + 1 is deduced from the location at time i. This idea yields 
an expected 0(n log n) complexity [BT93, pKS92 if the points are inserted in 
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a random order. The drawbacks of this approach are the following: the location 
structure consists of the history of the construction and thus strongly depends 
on the insertion order, and the additional memory needed cannot be controlled. 
(The expected memory is proved to be 0(n) and is experimentally about twice 
the size of the final triangulation.) 

Mulmuley [Mul91 proposed a location structure independent of the insertion 
order. The structure has 0(log n) levels, each level being a random sample of 
the level below. At each level, the Delaunay triangulation of the points is 
computed, and the overlapping triangles at different levels are linked to enable 
location of new points. This structure has the advantage of being independent 
of the order of insertion, of ensuring an 0(log 2 n) location time for any point, 
and of allowing deletions in an easier way than the Delaunay tree | DMT92 |. 
However, the additional memory is still important and the location structure is 
not especially simple. 

In 1996, Miicke, Saias and Zhu | MSZ96 | proposed a very simple structure 
to handle triangulation of random points. The structure reduces to a random 
subset of \fn points, and pointers from these points to an incident triangle in the 
Delaunay triangulation. A new point is located by finding the nearest neighbor 
in the sample by brute force, and walking in the triangulation. For evenly 
distributed points, the expected complexity of the algorithm is O(na) with a 
small constant, which makes it competitive with many 0[n log n) algorithms. 
But for some data (for example points on a parabola) the complexity increases 
to 0{ni). 

Overview 

Our approach uses a structure with levels similar to Mulmuley, but with 
simple relations between levels. This allows better control of the memory over- 
head. The transition between two levels is not direct as in Mulmuley, but uses 
a march similar Miicke, Saias and Zhu to locate point in triangulations. 

In Section ^| we present the algorithm, in Section |^ we prove that the ex- 
pected complexity of constructing the Delaunay triangulation is 0{n log n) . The 
parameters of the data structure are then tuned to minimize the constant in the 
case of random points and are shown to yield an excellent behavior in Section 
|], we pay special attention to the comparison with the method of Miicke, Saias 
and Zhu. Finally we give some implementation remarks and practical results in 
Section |[ 



2 Algorithm 

Let S be a set of n sites in the plane. The aim is to compute the Delaunay 
triangulation T>T $ of S and to maintain it efficiently under insertions and dele- 
tions. 
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2.1 The location structure 



The algorithm uses a data structure composed of different levels. Level i con- 
tains the Delaunay triangulation VTi of a set of sites Si . 

The sets Si forms a decrea sing sequence o f random subsets of S based on a 
Bernoulli sampling technique [MR95, Mul94| : 

5 = 5 D5 1 D5 2 3...D S fc _i D S k 

Prob(peS l+1 |j>€Si) = ie]0,l[. 

a 

The data structure is fairly simple: it contains the points of S and the 
triangles of all the triangulations T>Ti. A point p 6 S such that p 6 Si C . . . C 
iSo and p ^ <Si+i is said to be a vertex of level i and has a link to a Delaunay 
triangle of VTj incident to p for all j for < j < i. A triangle of VTi has links 
to its three neighbors in T>Ti and to its three vertices. The number k of levels 
is not fixed; for each point random trials decide its level, and the point with 
highest level determines k. 



2.2 Location of a query 

For the location of a query q, we start at a known vertex Vk+i of the highest 
level k. Then we search for Vk, the vertex of VT k nearest to q. Since Vk is also 
a vertex of T>Tk-i, we search for Vk-i, the nearest neighbor of q in VTk-i, 
starting at Vk ■ The search is continued descending the different levels. At each 
level i, the nearest vertex Vi of q in T>Ti is determined. 
At level i the search of Vi is carried out in three phases: 

• First phase: from Vi+±, we have a link to a triangle of T>Ti having u J+ i 
as vertex. All triangles incident to Vi+i are explored to find the triangle 
containing the segment Vi+iq. 

• Second phase: all the triangles of T>Ti intersected by Vi + \q are visited, 
walking along the segment Vi+iq up to the triangle ti that contains q. 

• Third phase: using neighborhood relationships between triangles, we will 
traverse few triangles of VTi from ti to find u,. If vv'v" are the three 
vertices of ti, and, without loss of generality, v is closer to q than v' and 
v" , then Vi is cither v or it lies in the disk of center q and passing through 
v (shaded on Figure |I]a); thus the search for Vi has to be done only in 
the direction of the neighbors of ti through the edges vv' and vv" and the 
neighbor through the edge v'v" can be ignored (the portion of the shaded 
disk in that direction is inside the disk through vv'v" which is empty) . 
For each such triangle, the distance to the new vertex is computed and 
the algorithm maintains the closest visited vertex. For a visited triangle 
ww'w" such that w is the nearest to q among ww'w" the neighbor triangle 
through edge ww' (resp ww") will be visited if angle qww' is smaller than 
f (Figure ||b). 

Figure He show the triangles visited by the different phases of the search. 
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Figure 1: Search for Vj . 



2.3 Updates 

Because of its simplicity, the data structure is fairly easy to update. Maintaining 
it dynamically provides a fully dynamic triangulation algorithm. The links 
between the different levels do not use any complicated data structure simply 
vertices know a triangle at all levels in which they appear. 

To delete a point from S, just delete the corresponding vertex at all the 
levels where it appears, which can be done in time sensitive to d the degree of 
that vertex. On average d = 6 and thus some of the following algorithms can 



be used. A complicated algorithm AGSS8E] of deterministic complexity 0(d), 
a simple randomized 0{d) algorithm |Chc86 can be used or simpler solutions 



of complexity 0(d\ogd) or even 0(n 2 ) may be good in practice. 

Inserting a point in S reduces to locating the new point at all levels, com- 
puting its level i and inserting the new vertex at all levels j, < j < i (which 
is sensitive to the degree of the new v ertex o nce the location is done). The 



insertion using the standard algorithm |Law77| 



3 Worst-case randomized analysis 

The analysis will rely on the randomization in the construction of the random 
subsets Si and the points of S are assumed to be inserted in a random order. 
In this section, no assumption applies to the data distribution, which can be in 
the worst case. As usual in theoretical computational geometry, we make only 
an asymptotic analysis and give rough upper bounds for the constants. In the 
next section, parameter a will be tuned to get a tight constant in the special 
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case of evenly-distributed points. 

Let S be a set of n points organized in the structure described in Section 
and q a point to be inserted in S. Since we have assumed a random insertion 
order, q is a random point of S U {q}. 

We denote m = |<Sf| and TZi — SiU {q}. 

Notice that, thanks to the random insertion order, TZi is a random subset of 
size rii + 1 of TZi-i and q is a random element of TZi. 

The cost of exploring all the triangles incident to Vi+i at the hrst phase of 
the march of level i is the degree of Vi+i in T>Ti. The cost of the second phase 
is the number of triangles intersected by segment Vi + \q. The cost of the third 
phase is the number of candidate vertices visited during the search of Vi from 
U. 

Lemma 1 The expected degree of Vi in T>T i-\ is 0(1). 

Proof Let J\fJ\f be the nearest neighbor graph of TZi: that is, the vertices 
of AW are the points of TZi, and q,v £ TZi define an edge of AW if and 
only if v is the nearest neighbor of q (denoted by v = NN(q)) or q is the 
nearest neighbor of v in TZi. AW is well known to be a subgraph of T>Tn i , 



the Delaunay triangulation of TZi, and to have maximum degree 6 [PY92 

We denote by d£^-._ (u) the degree of v in Z>7~j_i, and by E vE fc. c r q y the 
expectation when v is chosen uniformly in TZi C {q}. Then we have 

Evening} fer^ («)) = ^e^-icfg} (^r.^W) < 6 

notice that d^j-j-. (v) is a random variable; result holds since TZi and 
TZi-i C {q} are random subsets of TZi-% and that the average degree of 
a vertex in a triangulation is less than 6. 

But even if q is a random point in TZi, the vertex V{, the nearest neighbor 
of q in TZi, is not uniformly random. 



E, 



(fi Ti _ x (NN(q))) = E^J^d^r^NNiq))^ 

1 M \veKi qe{p;v=NN(p)} ) 
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Lemma 2 Given w ETZi, the expected number of vertices q of TZi such that w 
belongs to the disk of center q and passing through the nearest neighbor of q in 
IZi+i is less than 6a. 

Proof Let w £ IZi and let w = go, <?i , <?2 • • • Qk be the points of TZi lying 
in a section of angle ^ having apex w sorted by increasing distance to w. 
Clearly, a disk of center qi passing through qj (j < I) cannot contain w 
and thus, if q = qi, a necessary condition for w to be in the disk having as 
diameter the segment defined by q and the nearest neighbor of q in is 
that no point of {qo, ■ ■ ■ qi—i} is in the sample IZi+i which has probability 
(1 — 

Using six sections around w to cover the whole plane, and summing over 
the choice of q € TZi we get the claimed result. Notice that the disk of 
center q and passing through the nearest neighbor of q contain the disk of 
diameter the line segment defined by these two points, and thus the bound 
apply also to that circle. ■ 



Lemma 3 The expected number of edges ofDTi intersecting segment qvi+\ is 
0(a). 

Proof Let e be an edge of T>Ti intersecting segment If e does not 

exist in "DT^, it means that e is an internal edge of the region retrian- 
gulated when q is inserted in 2X7";. Since q is a random point in TZi, the 
expected number of such edges is 3 since it equals the average degree of q 
in TZi minus 3. 

If e exists in DTn i , one end-point w of e must belong to the disk of diameter 
qVi+i, denoted disk[(7« i+1 ], (otherwise any disk through the end-points of e 
must contain q or i> i+1 and e cannot belong to PT^J. 

The expected number of edges of VTn i intersecting disk[gu J+ i] is bounded 
by the sum of the degrees of the vertices in disk [qi)i+ 1] 

i?(#{e S having an end-point £ Kill disk 

1 <?eTC ; u^-Rindisk [qv i+ i] 

= "RTF E d VTK- H \ii e n *\ w e disk 

- Tul E d T>T- R . ( w ) Ga using Lemma| 

< 36a using the bound of 6 on the average degree of w 

Notice that Lemma ^] was established for a fixed w and a random q which 
allows to use it inside the sum over w. Thus we get a total expected cost 
for the march bounded by 36a + 3. ■ 
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Lemma 4 The expected number of triangles of T>T~i visited during the search 
for Vi from ti is 0{a). 

Proof All the triangles t examined in phase 3 have a vertex in the disk of 
center q passing through Thus we can argue similarly as in Lemma 

[| denoting disk \ c qv i+ i] the disk of center q through Vi + \: 

E(ff={t 6 VTfii having an end-point £ disk | c q«i + i]}) 

< iE E W«o 

^ TfT\ J2 d °vT^M \{q G ^> G 7^ n disk | c9 Wj + i]}| 

- TlTT E rf I?r TC .( w ) 6a using Lemma | 

< 36a using the bound on the average degree of to 



Theorem 5 TVie expected cost of inserting n point in the structure is 0(a log a n) 

Proof By linearity of expectation, Lemmas [l], || and ^ prove that the 
expected cost at one level is O(a). Since the expected height of the structure 
is log Q ro, we get the claimed result. (The analysis is similar to the ananlysis 



for skip lists |MR95|.) 



Theorem 6 The construction of the Delaunay triangulation of a set of n points 
is done in expected time 0(an\og a n) and 0{— zyn) space. The expectation is 
on the randomized sampling and the order of insertion, with no assumption on 
the point distribution. 

Proof Easy corollary of Theorem || ■ 



4 Tuning parameters 

We have proved that our structure is worst case optimal in the expected sense 
for any set of points. In this section, we will focus on more practical cases, and 
tune the algorithm to be optimal on random distribution. In that case, many 
events such as that a point has high degree and that it is the nearest neighbor 
of a random point can be considered as independent. 
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Figure 2: Number of orientation tests in phase 1 



4.1 Phase 1 

We can assume that, d^,^-. (fj+i) = 6 (and not only < 36 as proved in Lemma 
|l|). And thus if the turn around Vi+i is done in clockwise or counterclockwise 
direction depending on the position of segment Vi+\q with respect to the starting 
triangle, and assuming that this position is random around Vi+i the expected 
number of orientation tests is 3. Figure ^| shows the different cases to average, 
the edges Vi+iw such that an orientation test Vi+iwq is performed are indicated, 
for a typical degree 6 vertex in the triangulation. 



4.2 Phase 2 



Bose and Devroye [BD95] proved that the expected number of edges of a Dc- 
launay triangulation of random points crossed by a line segment of length I is 
0(ly/j) where 7 is the point density. Our experiments shows that the constant 
is 2. 

The expected number of points in the disk of center q passing through 
is a — 1. Indeed, if the points of IZi are sorted by increasing distance from q, 
v i+ i is the first point in lZ i+ i, thus the number of points in the disk is k with 
probability (1 — ^) fe ^, and the expected number is ^ XX^ — ^) fc = a — 1. Thus 
if I is the length of qv i+ i the density of points in X>T \ is 

Thus we conclude that the expected number of edges of TXT i intersecting 
segment qv l+1 is 21 ^f-^ = 2^2. 

For each edge ww' crossed, two orientation tests are performed: if w is the 
newly examined vertex, orientations of triangles wqvi + i and qww' are computed. 

We have to point out, that in the orientation tests of kind wqvi+i, the edge 
qvi + \ remains constant, and thus some computations do not need to be done 
for each test. 



4.3 Phase 3 



Phase 3 is more difficult to analyze precisely, but a rough bound is that the 
number of candidate vertices examined (with shortest distance) is less than two 
and that we examine less than 8 triangles in total. 

In fact, we modified phase 3, instead of really searching for Vi, the nearest 
neighbor of q in Si, we just define Vi as the nearest among the three vertices 
of tj. Thus this modified phase 3 reduced to three distance computations and 
two comparisons. 



4.4 Tuning a 

We will count more precisely the number of operation needed to evaluate our 
primitives. More exactly, we count the number of floating point operations 
(f.p.o.) without making diistinctions between additions, subtractions or multi- 
plications. 

The total evaluation at a given level is 3 + ^= orientation tests involving 

qvi+i, ^= other orientation tests and 3 distance computations. 

Orientation tests always using points q and Vi + i can be done using 5 f.p.o. 
to initialize plus 4 f.p.o. for each test. Other orientation tests need 7 f.p.o. each, 
and square distance computations need 5 f.p.o. each. 

Thus the total cost in terms of number of f.p.o. at level i is 

5 + 4(3 + + 7^= + 5 • 3 = 32 + 6.2 Va. 
V" V 71 " 

Since the number of level is log Q n — j°| 2 n a we get a cost of cq (n) = (32 + 

j°| 2 n a which is close to its minimum ( S [13.3 log 2 n, 141og 2 n]) for 
a € [18,90], with the minimum occuring for a ~ 40. 



4.5 Comparison with [ |M$Z9q 



Similar counting of f.p.o. in Miicke et al. algorithm, using a random sample of 
{3\/n points, produces a cost of 

n n 
0n3 . _ /3n3 _ „ »i— ^ _ ,1—1 6.2 



CMSz{n) = 5 + 4(3 + ^) + 7^ + 5/3^ = 17 + jfc ^= + 5/3 

V 71 " V 71 " Vvp . 

which is close to its minimal value for 0.5 < p < 1. 

As shown by the comparison of the two curves in Figure [| our method 



is potentially much better than |MSZ96|, even for a small number of points 



However, this method to analyze our approach hides the discontinuity of the 
cost, since the effective number of levels is necessarily an integer. To have a 
better comprehension of what happens for a small number of points, we can 
draw the cost of inserting a point in a structure having a fixed number of levels. 
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number of points 



Figure 3: Comparison of number of floating point operations between cq(ji) and 
CMSz{n) for a = 40 and j3 = 1. 
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The classical walk from a random point in the structure costs 



c wa ik{n) = 5 + 4(3 + ^) + 7^ = 17 + 6.2^ 

/7T 



which is also the cost of inserting in our structure up to the time a second level 
is created. 

When fc levels have been created, the cost is 

c k (n) = c waik {^k) + 15fc + fc • c walk (a) 

We can alternatively mix this multilevel approach with Miicke et al's. sam- 
pling at the first level of the structure. In that case, the cost is 

c* k (n) = c AI sz {^) + 15fc + fc • c wa i k (a) 



This comparison (see Figure [|) shows that [MSZ96 (c\{n)) becomes better 
than the simple march (ci(n)) for n > 40. The two level structure (c2(ro)) 
becomes better than the single level structure (ci(n)) for n > 180 and better 
than [MSZ96| (c^(n)) for n > 600. The main information is that the structure 



presented in that paper should be significantly better than [MSZ96| for 10000 < 



5 Implementation 

5.1 Deletion 

The above structure supports insertions and queries as explained above, but also 
deletions. Since there is no complicated data structure to maintain, deletions 
can be handled by just deleting the removed poi nt at ea c h level w here it appears. 
This can be done in output-sensitive time Che87 , AGSS89|] , and thus the 



deletion of a random point is done in expected constant time since a point 
appears at an expected constant number of levels and its expected degree fc is 
also constant. 

From a practical point of view, and to keep the simplicity of the algorithm, 
a simpler suboptimal algorithm should be preferred. It can be done in 0(fc 2 ) 
time, for example by flipping to reduce the degree of the deleted vertex to 3, 
and flipping again to restore the Delaunay property. Another simple algorithm 
consists in finding the Delaunay triangle incident to an edge of the hole in O(fc) 
time which also yields an 0(fc 2 ) time algorithm. Both algorithms are efficient in 
practice and needs only few micro-seconds (about 30 in a random triangulation) 
to delete a point once it had been localized. 

5.2 Arithmetic degree 

The algorithm above is designed to make a parsimonious use of high degree tests 



[TLP96|. More precisely, the location phase uses only orientation tests on three 
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Figure 4: Comparison of number of floating point operations between Cfc(n) and 
c\(n) for a = 40. 
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points in phases 1 and 2, and distance computation and angle comparisons with 
| in phase 3. All these tests are degree 2 tests. Clearly, updates need to use 
in-circle tests which are of degree 4. 

An alternative to phase 3 should have to use in-circle tests to limit the 
explored triangles in VT \ to those whose circumcircle contains q. Such variant 
may explore fewer triangles and be easier to analyze, but may use more degree 
4 tests. 

5.3 Robustness issues and degeneracies 

Degeneracies are solved by handling special cases: if two points have the same 
coordinates, then the insertion is not done, if four points are cocircular, then 
the last point inserted is considered as inside the disk defined by the others. 

We use exact arithmetic for 24 bits integers, and thus coordinates of our 
points are integers in range [—16777216, 16777216] (up to a multiplication by a 
power of 2). Using this restricted kind of data, double precision computation is 
exact on degree 2 tests and almost never leads to precision problems on degree 
4 predicates. Nevertheless, the exactness of all computations are verified by an 
arithmetic filter and exact computation is performed if needed. 

5.4 Code parameters 

The following parameters can be specified: 

• maximal number of levels 

• a the ratio between two levels 

• the minimal number of points to use the higher level for point location 

• the minimal number of points to use MSZ sampling at one of the higher 
levels 

• (3 the constant for the size of MSZ sample. 
Our default parameters are 

• number of levels unlimited 

• a = 30. 

• minimal size to use hierarchy is 20. 

• minimal size to use MSZ is 20. 
. (3 = 1. 

We found that the code is relatively insensitive to the parameters. For reason- 
able changes of these parameters, (up to a factor 2) the computation time is not 
greatly affected. Using these configuration parameters, our code can be used to 
run 
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the usual walk algorithm (only one level and minimal size for MSZ=oo), 
the Miicke et al. algorithm [ MSZ96[ (only one level), 



the hierarchical algorithm described in this paper (minimal size for MSZ=oo), 



the mixed method suggested in Section 4.5 (default parameters above) 



5.5 Experimental results 

5.5.1 Data sets 

We claim that our algorithm performs well on random point sets, and has ac- 
ceptable worse case complexity. To illustrate this fact, we will test it with the 
realistic and degenerate data sets. For each kind of data, we used sets of size 
5,000, 50,000 and 500,000 points. The coordinates are random on 24 bits and 
the constraints such that the points are on a parabola are verified, up to the 
rounding arithmetic errors. 

• random: points evenly distributed in a square. 

• ellipse: points evenly distributed on an ellipse. 

• ellipse2: 95% points evenly distributed on an ellipse plus 5% points evenly 
distributed in a square. 

• circle: points evenly distributed on a circle. 

• parabola: points evenly distributed on a parabola, 

If the circle and parabola examples can be considered as pathological inputs, 
the ellipse and ellipse2 examples are more realistic, Delaunay triangulation of 
points distributed on a curve occurs in practical applications, for example in 
shape reconstruction (see Figure ^|) . 

5.5.2 Results 

Following results are obtained on a Sun-Ultral 200 MHz. The code is written 
in CH — h and compiled with AT-T compiler with optimizing options. Time has 
been obtained with the clock command and is given in seconds. The time which 
is measured is just the Delaunay computation; it does not take into account the 
time for input or output. 

Figure [(] gives the computation times for execution of the code with the 



different parameters described in Section 5.4. Since it is the same code, the 
low level primitives such as in-circle tests or the walk in the triangulation are 
identical and it provides a fair comparison between the different methods. 

The last column is always the fastest method. It is significantly better 
than MSZ for very large sets of random points, and the difference is even more 
important on data set ellipse2 which is representative of real applications. 
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parabola ellipse ellipse2 



Figure 5: Data sets. 
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Figure 6: Running times 



5.5.3 Comparison with other software 

We have compared with some Delaunay softwares available on the WWW: 

• qhull by Bradford Barber and Hannu Huhdanpaa, duality with 3D convex 

hull |BDH93|] (available at 

jxttp : / / www, gcom . umn . cdu /locate /qhull ) . 



• div-conquer by Jonathan Shcwchuk, divide and conquer [She96 

• sweep by Jonathan Shewchuk, plane sweep 

• incremental by Jonathan Shewchuk, incremental with Miicke et al. lo- 
calization. 

These three codes supports exact arithmetic on double (available at 
frittp://www. cs.cmu.edu/~quake/triangle.rcscarch.htmj ). 



Dtree Delau nay tree structurc| BT93 | (time includes input) 
(available at http:/ /www.inria.fr/prisme/logiciel/del-tree.html). 



• hierarchy this paper, mixed with MSZ. 

The execution times in seconds are in Figure ^. Our method is significantly 
faster than the other incremental method, especially in the ellipse cases. Our 
method is about 50% slower than the divide and conquer algorithm. 
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Figure 7: Comparisons with other softwares 



6 Conclusion 

We proposed a new hierarchical data structure to compute the Delaunay trian- 
gulation of a set of points in the plane. It combines good worst case randomized 
complexity, fast behavior on real data, small memory occupation and dynamic 
updates (insertion and deletion of points). 

Referring to Su and Drysdale []SD97 1 st udy of several techniques and our 



comparisons with Shewchuk implementation [She96| of some of these techniques, 
we have shown that our implementation is competitive with other approaches on 
random data. Furthermore, we can prove that the performances remains good 
on pathological inputs. Finally, one of the main advantage of this algorithm is 
to allow a dynamic setting. 

The main idea of our structure is to perform point location using several 
levels. The lowest level just consists of the triangulation, then each level contains 
the triangulation of a small sample of the levels below. Point location is done 
by marching in a triangulation to determine the nearest neighbor of the query 
at that level, then the march restart from that neighbor at the level below. 



Location at highest level is done using [MSZ96] which is efficient for small set 
of points. 

One characteristics of the structure is that best time performance is obtained 
with a ratio of about three per cent between two levels, which yields to few levels 
(three or four typically) and a small memory occupation. The structure is simple 
and does not need additional features such as buckets. 

Such structure can be generalized to other problems. The two main ingre- 
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dients of the proofs are bounds on the maximal degree of the nearest neighbor 
graph and the expected degree of a random vertex in the Delaunay triangula- 
tion. The first generalizes well in higher dimension, while the second becomes 
an data sensitive parameter (constant for random points, n^ d_1 " 2 1 in the worst 
case). A generalization for computing the trapezoidal map can also be done. 

Code 

A demo version compiled for Sun Solaris and SGI is available at 
http: / /www. inria.fr/prismc/logicicls/del-hierarchy/. 



Acknowledgement 

The author would like to thank Herv Brnnimann, Jonathan Shewchuk, Jack 
Snoeyink and Mariette Yvinec for helpful discussions and careful reading of this 
paper. 



References 



[AGSS89] A. Aggarwal, Leonidas J. Guibas, J. Saxe, and P. W. Shor. A linear-time 
algorithm for computing the Voronoi diagram of a convex polygon. Discrete 
Comput. Geom., 4(6):591-604, 1989. 

[BD95] P. Bose and L. Devroye. Intersections with random geometric objects. 

Technical report, School of Computer Science, McGill University, 1995. 
Manuscript. 

[BDH93] C. B. Barber, D. P. Dobkin, and H. Huhdanpaa. The Quickhull algorithm 
for convex hull. Technical Report GCG53, Geometry Center, Univ. of Min- 
nesota, July 1993. 

[BT86] Jean-Daniel Boissonnat and Monique Teillaud. A hierarchical representa- 
tion of objects: The Delaunay tree. In Proc. 2nd Annu. ACM Sympos. 
Comput. Geom., pages 260-268, 1986. 

[BT93] Jean-Daniel Boissonnat and Monique Teillaud. On the randomized con- 
struction of the Delaunay tree. Theoret. Comput. Sci., 112:339-354, 1993. 

[Che86] L. P. Chew. Building Voronoi diagrams for convex polygons in linear ex- 
pected time. Technical Report PCS-TR90-147, Dept. Math. Comput. Sci., 
Dartmouth College, Hanover, NH, 1986. 

[Che87] L. P. Chew. Constrained Delaunay triangulations. In Proc. 3rd Annu. ACM 
Sympos. Comput. Geom., pages 215-222, 1987. 

[DMT92] Olivier Devillers, Stefan Meiser, and Monique Teillaud. Fully dynamic De- 
launay triangulation in logarithmic expectedtime per operation. Comput. 
Geom. Theory AppL, 2(2):55-80, 1992. 

[GKS92] Leonidas J. Guibas, D. E. Knuth, and Micha Sharir. Randomized incremen- 
tal construction of Delaunay and Voronoi diagrams. Algorithmica, 7:381- 
413, 1992. 



18 



[Law77] C. L. Lawson. Software for C 1 surface interpolation. In J. R. Rice, editor, 
Math. Software III, pages 161-194. Academic Press, New York, NY, 1977. 

[MR95] R. Motwani and P. Raghavan. Randomized Algorithms. Cambridge Univer- 
sity Press, New York, NY, 1995. 

[MSZ96] Ernst P. Miicke, Isaac Saias, and Binhai Zhu. Fast randomized point lo- 
cation without preprocessing in two- and three-dimensional Delaunay tri- 
angulations. In Proc. 12th Annu. ACM Sympos. Comput. Geom., pages 
274-283, 1996. 

[Mul91] K. Mulmuley. Randomized multidimensional search trees: Dynamic sam- 
pling. In Proc. 7th Annu. ACM Sympos. Comput. Geom., pages 121-131, 
1991. 

[Mul94] K. Mulmuley. Computational Geometry: An Introduction Through Ran- 
domized Algorithms. Prentice Hall, Englewood Cliffs, NJ, 1994. 

[PY92] M. S. Paterson and F. F. Yao. On nearest-neighbor graphs. In Proc. 19th 
Internat. Colloq. Automata Lang. Program., volume 623 of Lecture Notes 
Comput. Sci., pages 416-426. Springer- Verlag, 1992. 

[SD97] P. Su and R. Drysdale. A comparison of sequential Delaunay triangulation 
algorithms. Comput. Geom. Theory Appl, 7:361-386, 1997. 

[She96] J. R. Shewchuk. Triangle: engineering a 2d quality mesh generator and 
Delaunay triangulator. In First Workshop on Applied Computational Ge- 
ometry. Association for Computing Machinery, May 1996. 

[TLP96] R. Tamassia, G. Liotta, and F. P. Preparata. Robust proximity queries 
in implicit Voronoi diagrams. In Proc. 8th Canad. Conf. Comput. Geom., 
page 1, 1996. 



19 



